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TRAJECTORY DESIGN FOR A CISLUNAR CUBESAT 
LEVERAGING DYNAMICAL SYSTEMS TECHNIQUES: THE 
LUNAR ICECUBE MISSION 


Natasha Bosanac* Andrew Cox! Kathleen C. Howell: and David C. Folta: 


Lunar IceCube is a 6U CubeSat that is designed to detect and observe lunar volatiles from a 
highly inclined orbit. This spacecraft, equipped with a low-thrust engine, will be deployed 
from the upcoming Exploration Mission-1 vehicle in late 2018. However, significant un- 
certainty in the deployment conditions for secondary payloads impacts both the availability 
and geometry of transfers that deliver the spacecraft to the lunar vicinity. A framework that 
leverages dynamical systems techniques is applied to a recently updated set of deployment 
conditions and spacecraft parameter values for the Lunar IceCube mission, demonstrating 
the capability for rapid trajectory design. 


INTRODUCTION 


With the emergence and increased development of miniaturized satellite technologies, CubeSats offer an 
alternative platform for unmanned exploration of cislunar space and, eventually, the solar system.!:*»> In 
fact, their reduced cost and development time in contrast to larger, traditional spacecraft supports a wide 
array of organizations conducting scientific observation, technology demonstration and space exploration 
missions. Continued advancement in these small satellite technologies demands the innovative design of 
trajectories for spacecraft that require fewer resources, achieve more complex mission goals and visit even 
farther destinations. Furthermore, small spacecraft typically ride as secondary payloads onboard a launch 
vehicle, producing significant uncertainty and variability in the launch date and deployment status. However, 
such uncertainty in the deployment state of a small spacecraft can severely constrain and complicate the 
trajectory design process during both the development and operational phases of the mission. Accordingly, 
the construction of a framework that leverages dynamical systems techniques and enables efficient, rapid and 
well-informed trajectory design is warranted.*” 


To demonstrate the value of a trajectory design framework, consider the Lunar IceCube mission led 
by Morehead State University and supported by NASA Goddard Space Flight Center (GSFC), Busek and 
Catholic University of America. The mission concept is based on the use of a 6U CubeSat which is designed 
to detect and observe water near the Moon’s equatorial regions.” To achieve these scientific goals, an in- 
clined, low altitude lunar elliptical orbit is required during the observational phase of the mission. Currently, 
Lunar IceCube is scheduled for launch in late 2018 as a secondary payload onboard the upcoming Explo- 
ration Mission-1 (EM-1) vehicle. Following deployment close to the Earth, the Lunar IceCube spacecraft is 
deployed to a high energy trajectory that would naturally depart the Earth-Moon system. However, the space- 
craft features a low-thrust propulsion system that is currently estimated to deliver up to 1.15 mN of thrust 
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with an J, of 2500s. Finite duration maneuvers can be leveraged, in combination with the natural dynamical 
structures associated with the gravity of the Sun, Earth and Moon, to design a complex path that delivers the 
Lunar IceCube spacecraft to its final lunar science orbit.* However, the constrained departure conditions, 1.e., 
the launch date and deployment state, as well as the reduced propulsive capability, may limit both the geom- 
etry and persistence of transfers within the phase space. Furthermore, as these departure conditions evolve 
during both the development and operational phases of the primary mission, a complete trajectory redesign 
may be required, necessitating the development of a framework for rapidly constructing an initial guess for a 
transfer that may satisfy the mission and hardware requirements. 


To construct a trajectory design framework for a CubeSat with limited propulsive capability and con- 
strained departure conditions, techniques from dynamical systems theory are leveraged. In fact, natural 
dynamical structures, along with Poincaré mapping strategies, are employed to design a complex path that 
delivers the spacecraft from a highly energetic deployment state to the final lunar science orbit.* With applica- 
tion to the Lunar IceCube mission, analysis of existing point solutions reveals that the transfer trajectory can 
be split into three components: an Earth outbound segment, a phasing and energy adjustment segment that 
exploits solar gravity and, finally, lunar approach and capture.*’!° The first and final segments rely heavily 
on the low-thrust propulsion system to modify the spacecraft energy. However, given the limitations on the 
propulsive capability of the Lunar IceCube spacecraft, predominantly natural arcs are sought to link these two 
bounding segments. In addition, the phasing and energy adjustment segment exhibits the largest variability in 
terms of the available geometries and trajectory characteristics. To effectively design this transfer segment, a 
circular restricted three-body model of the Sun-Earth system is employed to supply insight into the geometry 
of natural dynamical structures that can be incorporated into the design of an end-to-end path. Furthermore, 
approximate bounds on the motion are established and transfer geometries are understood via analysis of the 
manifolds of libration point orbits.'!? These natural motions are visualized and linked to the deployment and 
final lunar science orbit bounding conditions via Poincaré mapping strategies, thereby enabling rapid and ef- 
ficient construction of a discontinuous initial guess. Depending on the deployment conditions, that frequently 
evolve throughout further development of the mission concept, the exact definition of the constructed maps 
may require updates; however, a general framework that leverages mapping techniques remains useful. A 
differential corrections process is then formulated in a low-thrust-enabled ephemeris model to identify an 
end-to-end trajectory that generally retains the characteristics of the initial guess.'° The resulting continu- 
ous solution can then be input to an operational-level model for further corrections or even an optimization 
algorithm. The constructed trajectory design framework, employed to link the highly energetic deployment 
state of the Lunar IceCube spacecraft to the final lunar science orbit, supports the identification of feasible 
transfer regions and rapid redesign. Implementation of this strategy eliminates the challenges associated 
with a grid search in a nonlinear and chaotic system. Moreover, this framework may be applicable to future 
CubeSat missions that must meet alternative mission goals, such as re-encountering the Moon with a specific 
energy and/or flight duration, attaining specific Sun-Earth or Earth-Moon orbits, or achieving a heliocentric 
trajectory that encounters an asteroid, as well as mission extensions. 


LUNAR ICECUBE MISSION OVERVIEW 


The Lunar IceCube mission leverages a 6U CubeSat developed by Morehead State University and sup- 
ported by NASA Goddard Space Flight Center (GSFC), Busek and Catholic University of America. Specif- 
ically, the primary objective for this mission is to observe the equatorial regions of the Moon, furthering 
knowledge of the location and transportation physics associated with water in various forms.” The scien- 
tific instruments onboard the Lunar IceCube spacecraft require that these measurements of lunar water and 
volatiles be collected from a low lunar science orbit. Currently, the nominal science orbit is highly elliptical 
with a perilune located over the equator at an altitude between 100 and 105 km. Furthermore, the orbit is also 
constrained to possess an inclination of approximately 90 degrees. The Lunar IceCube spacecraft, depicted 
in Figure 1 and currently estimated to possess an initial mass of 14 kg, will ride onboard the Exploration 
Mission-1 (EM-1) vehicle as a secondary payload. The spacecraft is also equipped with an iodine-fuelled 
low-thrust engine, a Busek Ion Thruster 3-cm (BIT-3) system, which is currently designed to deliver up to 
1.15 mN of thrust with an /,,, of 2500s and a propellant mass of approximately 1.5 kg.° Previously-developed 


Figure 1: Lunar IceCube preliminary spacecraft design. 


trajectory point designs that leverage this low-thrust engine have revealed distinctly different trajectories in 
a large design space that is difficult to examine through the construction of single point solutions in an 
ephemeris environment, particularly as the spacecraft hardware and mission parameters evolve.* Thus, a 
framework employing dynamical systems techniques is used to construct similar solutions based on insight 
into the natural structures within the gravitational environment of the Earth, Sun and Moon. 


DYNAMICAL MODELS 
Circular Restricted Three-Body Problem 


Motion within the Sun-Earth system is rapidly and reasonably approximated using the autonomous dy- 
namics of the Circular Restricted Three-Body Problem (CR3BP). In this dynamical environment, the motion 
of a massless spacecraft is modeled under the influence of the point-mass gravitational attractions of two 
primaries: the Sun and the Earth. To enable clear visualization and identification of particular solutions, the 
motion of the spacecraft is described using a rotating coordinate system. This frame, (x, y, z), rotates with 
the primaries as they encircle their mutual barycenter. In addition, position and velocity states locating the 
spacecraft are nondimensionalized. By convention, both the normalized distance between the Sun and the 
Earth and the mean motion of the primaries are unity. Mass quantities are nondimensionalized such that the 
masses of the Sun and the Earth are equal to 1 — py and yp, respectively. Using these nondimensionalized 
quantities, the equations of motion of the spacecraft, located at (x, y, z) in the rotating frame, are compactly 
written as: 
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(cx —1+ ys)? + y? + z?. When formulated in the rotating fants this eee Saami environment 
admits a constant energy-like integral labeled the Jacobi constant, Cy; = 2U — x? — y? — 27. Ata fixed 
value of this integral, an infinite number of trajectories exist within the Sun-Earth system. Any state along 
a time-varying solution is defined as prograde if the corresponding angular momentum vector at that instant 
possesses a +2 component, i.e., the spacecraft is traveling in a counterclockwise direction as viewed in 
the rotating frame of the CR3BP. Correspondingly, a state along a path that possesses a —z component of 
the angular momentum vector is labeled retrograde. Regardless of the direction of motion, these particular 
solutions exhibit one of four types of behavior - equilibrium points, periodic orbits, quasi-periodic orbits 
or chaotic motion - that may be leveraged to design a complex trajectory within the nonlinear gravitational 
environment of the Earth, Moon and Sun. 


Low-Thrust-Enabled Motion in an Ephemeris Model 


To support recovery of a representative end-to-end trajectory, equations of motion are formulated for a 
spacecraft traveling under the point mass gravitational influences of the Earth, Moon and Sun, as well as 
an additional acceleration due to a low-thrust propulsion system, leveraging a spacecraft-centered velocity- 
normal-conormal coordinate frame to describe the direction of the thrust vector.! First, consider N, bodies 
with dimensional mass /;, each assumed to be spherically symmetric and, therefore, modeled as a point 
mass. These bodies are located in an inertial reference frame, XY Z, with a fixed origin at point O. Then, 


the dimensional position vector locating the body of interest P3, at the inertial coordinates R3 = (X V7 y 
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with respect to the body P;, at the coordinates R; = (X;,Y;,Z;), is found as Rjz = (X — X;)X + (Y - 
Y;)Y + (Z — Z;)Z. To formulate the differential equations governing natural motion within the point mass 
ephemeris model, the acceleration of the body of interest, P3, is expressed relative to a celestial body P, 
in P,-centered J2000 inertial coordinates.’ When P3 represents a spacecraft, M3 < M,, M ;- Thus, the 
spacecraft is assumed to only negligibly impact the paths of the celestial bodies. Accordingly, additional 
equations of motion are not required for the VV. — 1 celestial bodies if ephemerides are available. 


Prior to augmenting the point mass ephemeris model equations of motion with the additional acceleration 
contributed by a propulsion system, the characteristics of a low-thrust engine are summarized. Specifically, 
a constant specific impulse low-thrust engine can be described by its thrust magnitude, 77;, specific impulse, 
Is», and the thrust direction, w. For this investigation, the iodine-fuelled BIT-3 system onboard the Lunar 
IceCube spacecraft is assumed to operate with the following parameters: T;, = 0.9 mN and J, = 2500s.° It 
is also assumed that the mass of the spacecraft, m(t), is initially equal to 14 kg. When the low-thrust engine 
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is activated, the spacecraft mass decrements by the following mass flow rate, ane — —5e 
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constant value for the available engine power is calculated as P = Att“spdo for a gravitational acceleration, 
go = 9.81m/s?.° This constant mass flow rate impacts the maximum time for operation of the low-thrust 


engine which is constrained by the initial propellant mass, assumed to be equal to 1.5 kg.* 


where the 


In this analysis, the direction of the thrust vector is described using a velocity-normal-conormal (VNC) 
coordinate system, which is often employed for spacecraft applications.” Specifically, the velocity direction, 
denoted by the V unit vector is parallel to the inertial velocity vector of the spacecraft relative to a reference 
body. Then, the conormal direction, C, is defined normal to the path of the spacecraft, while lying within the 
orbital plane. Finally, the normal unit vector, N,is perpendicular to the orbital plane and completes the right- 
handed coordinate system. The engine thrust direction, u, can then be projected into the VNC coordinate 
system and expressed using these three unit vector directions such that u = uyV tunN +ucC." The 
thrust vector, Wy wc, expressed in the VNC frame can subsequently be converted into the P,-centered J2000 
inertial coordinates to produce, Uxyz. This transformation is leveraged to incorporate a thrust vector that 
possesses a fixed direction in the VNC frame into the equations of motion for a point mass ephemeris model. 


The additional acceleration for a low-thrust engine is incorporated into the differential equations govern- 
ing the motion of a spacecraft in a point mass ephemeris model including the Sun, Earth and Moon. First, 
to mitigate ill-conditioning between the position, velocity and time quantities in the equations of motion, 
nondimensionalization is employed. Specifically, the characteristic quantities m*, t*, and /*, computed in the 
CR3BP for a primary system consisting of the bodies P, and P,, may be employed. These nondimensional 
quantities are denoted without a tilde. Then, the acceleration due to the low-thrust engine is expressed in 
vector form as 7;;/m/(t)t where the dimensional mass of the spacecraft, m, is employed, and the symbol 

7, represents the thrust magnitude nondimensionalized only in the length and time units for dimensional 
consistency such that Ty) = Ty, x (t*)*/l*. Then, the coefficient T;*/m is nondimensional. Once @ is 
transformed into the same coordinate frame as the acceleration of the spacecraft relative to a body Py, the 
additional acceleration term is straightforwardly added to the equations of motion. When modeling the mo- 
tion of a spacecraft relative to the Moon within a Sun, Earth, and Moon point mass ephemeris model, the 
low-thrust-enabled equations of motion are written as: 
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where the subscript FE identifies the Earth, s/c identifies the spacecraft, LZ corresponds to the Moon and S$ 
indicates the Sun. Since these three equations of motion feature the time-dependent mass of the spacecraft, the 


mass flow rate equation is also required and must be simultaneously integrated along with the spacecraft state. 
ge these differential equations are used to numerically integrate the augmented state of the spacecraft, 
Ke = |G, X.V7 ,m]|, where the nondimensional position and velocity components are expressed 
hive to the Moon in the Moon J2000 inertial coordinate system. At each epoch during the integration 
process, the position coordinates corresponding to each body, relative to a designated basepoint located at the 
Moon, are accessed using the Jet Propulsion Laboratory DE421 ephemerides via the SPICE toolkit.!?°!3!4 
Relative position coordinates are then used to integrate the motion of a spacecraft relative to the Moon in 
the low-thrust-enabled four-body point mass ephemeris model. Alternatively, for motion predominantly in 
the Earth vicinity, the equations of motion can be written relative to the Earth for better performance during 
numerical integration. 


DIFFERENTIAL CORRECTIONS TO RECOVER AN END-TO-END TRAJECTORY 


To construct a continuous trajectory for a low-thrust-enabled spacecraft traveling within a point mass 
ephemeris model of the Earth, Moon and Sun, a variable-time multiple shooting algorithm is employed.!° 
The general procedure for a multiple shooting algorithm involves intermediate nodes, described by a set of 
free variables, that are simultaneously integrated to produce a sequence of arcs and, then, iteratively corrected 
to recover a continuous trajectory that satisfies a set of constraints. To demonstrate the formulation of the 
multiple shooting algorithm employed in this application, the conceptual representation in Figure 2 serves as 
a reference. In this illustration, the path of a spacecraft relative to a central body is sought as it travels under 
the gravitational influence of the central body and any additional perturbing bodies, as well as the action of a 
low-thrust engine. 


Consider n, nodes, indicated via blue circles in Figure 2, to produce a sequence of arcs that, when inte- 
grated forward in time, form the initial guess for a trajectory. The 2-th node is completely specified by the 
vector 5; = [X;, Yi, Zi, x Y.. Zi, ie | eden: ieee In particular, the description of each node includes: 
its nondimensional position, (X;, Y;, Z;), and velocity, (Xi, Yj, 7. relative to body P, in a P,-centered 
J2000 inertial frame; the spacecraft mass, m,;, to ensure mass continuity across arcs that employ low-thrust; 
the epoch at each node, Tirode,i, expressed relative to a fixed initial epoch and nondimensionalized via the 
characteristic time quantity, “, of an appropriate P, — P, system; and the nondimensional integration time, 
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Figure 2: Illustration of variable time multiple shooting formulation to target a continuous trajectory for a 
low-thrust-enabled spacecraft in an ephemeris model. 


Tint,i- Together, these variables completely represent each node. However, to completely describe the arc as- 
sociated with a given node, a constant binary parameter is used to distinguish between a natural arc, displayed 
in blue in Figure 2, and a low-thrust arc, indicated in brown. Although this parameter must be predefined, 
the integration time along an arc is allowed to vary towards zero, thereby enabling a low-thrust arc to be 
essentially eliminated, if necessary, to recover a continuous path. For a low-thrust-enabled arc, the three 
components of the unit thrust direction vector, wy, uy, and uc, must be specified. To identify arcs over 
which the low-thrust engine is activated, as well as the basepoint body used to express the position and ve- 
locity vectors, constant vectors with integers encoding this information are employed. Using this description, 
each of the arcs along a trajectory is completely defined. 


Multiple arcs, computed by integrating a set of nodes forward in time, are iteratively refined until they sat- 
isfy a sequence of constraints that produce a continuous trajectory that leverages both natural and low-thrust- 
enabled motion. First, connectivity between two neighboring arcs requires state, mass and time continuity. 
In particular, the position, velocity and mass at the end of the 2-th arc, integrated forward for a time Tjnz,i, 
must equal the position, velocity and mass at the 2 + 1-th node to within an acceptable tolerance. In addition, 
since the ephemerides of the central and perturbing bodies are time-dependent, the epoch, Thode, + Tint,is 
at the end of the z-th arc must equal the epoch, Tiode,it1, at the 2 + 1-th node; this condition is labeled a 
time continuity constraint.’ Furthermore, to ensure that the integration time along the i-th arc is positive, an 
inequality constraint is converted to an equality constraint that leverages a slack variable, 3;. This constraint 
is expressed in the form Tint; — Be = 0. Then, for any low-thrust-enabled arcs, an additional constraint is 
required: the thrust vector components, uy, un, and uc, must produce a unit vector. Thus, the unit thrust 
vector constraint is written as u?, + u%, + uz, — 1 = 0. The required constraints for the i-th segment are then 
summarized as follows: 
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and over the segments for which the low-thrust engine is activated, uj, , + uy, + ug, — 1 = 0. Together, 
these constraints are simultaneously employed at each step of an iterative corrections algorithm to recover 
a continuous trajectory for a low-thrust-enabled spacecraft under the influence of a point mass ephemeris 
model of the Earth, Moon and Sun.!° In practice, additional constraints may need to be enforced at selected 
nodes along a trajectory. Such constraints, which may link to mission requirements or may be required for 
analysis, involve the epoch, spacecraft mass, inertial state components, apse conditions, energy via the Jacobi 
constant, rotating frame state components or orbital elements. When required, such constraints are added to 
the continuity, positive time and unit thrust vector constraints in the multiple shooting algorithm implemented 
within MATLAB, and their derivatives may be computed using any combination of analytical expressions and 
numerical approximations. !° 


TRAJECTORY DESIGN FRAMEWORK 


By leveraging known natural dynamical structures that exist in a simplified, autonomous dynamical model, 
a rapid and well-informed process for constructing a trajectory for a low-thrust-enabled small satellite is 
demonstrated. In particular, the application of dynamical systems techniques enables the identification of 
known solutions that exist within a simplified natural gravitational environment. This insight offers a scheme 
to predict the bounds on the motion, while also supporting the construction of an initial guess for the trajec- 
tory prior to corrections in a operational-level modeling environment.* The capability for a rapid and guided 
exploration of the trajectory design space is particularly beneficial for CubeSat missions that possess limited 
propulsive capability and, as a secondary payload, are subject to significant uncertainty in their deployment 
conditions. Furthermore, the construction of a framework for CubeSat trajectory design, as demonstrated 
for the Lunar [ceCube mission, enables rapid redesign as the deployment and spacecraft hardware param- 
eters evolve during both the development and operational phases of the mission. This framework involves 
decomposing the trajectory into three segments: 


e An Earth outbound segment, originating at deployment 
e A phasing and energy adjustment segment that exploits solar gravity, and 


e A lunar approach and capture portion. 


To illustrate this trajectory sequence, consider a sample point solution, originally designed by David Folta 
at NASA GSFC for a thrust magnitude of 1.2 mN and an earlier set of deployment conditions (from August 
2015), as depicted in the top plot within Figure 3.*’!! In this figure, blue portions along the trajectory indicate 
natural motion while green curves correspond to arcs where the low-thrust engine is activated. In the three 
plots at the bottom of Figure 3, the arcs comprising each of the three segments are highlighted, while portions 
of this sample solution that are not included in the segment of interest are colored gray. Of these three 
segments along the trajectory sequence, the focus of this investigation is the phasing and energy adjustment 
arc, which exhibits the largest variability in terms of the available geometries and path characteristics. To 
explore and identify natural dynamical structures that may offer a pathway to link a high energy deployment 
state to a lunar capture arc, a Sun-Earth CR3BP model is leveraged, along with Poincaré mapping strategies. !° 
This framework is employed to demonstrate the construction of an initial guess for a new set of deployment 
conditions via dynamical systems techniques and insight. This initial guess trajectory is then corrected in a 
point mass ephemeris model that incorporates the additional acceleration of a low-thrust engine. The resulting 
trajectory supplies an initial guess for corrections and optimization within a higher fidelity dynamical model 
and, eventually, an operational-level modeling environment. 
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Figure 3: Illustration of the trajectory sequence employed in this investigation, demonstrated using a point 
solution for the original deployment conditions (top). Portions of the trajectory comprising each of three 
segments are isolated in the inset images (bottom). 


Earth Outbound Segment 


Given the original trajectory design, as illustrated in Figure 3, adapting to a new deployment state is 
nontrivial. For example, new deployment conditions are introduced and are current as of December 2016. 
Simulation of these most recently updated deployment conditions reveals that, naturally, the associated tra- 
jectory quickly departs the Earth vicinity. For a deployment epoch date of October 7, 2018, the initial state 
corresponds to a highly energetic path when integrated forward in time in a point mass ephemeris model and 
incorporating a lunar encounter. In fact, the associated path is depicted in nondimensional coordinates in 
Figure 4 in (a) the Sun-Earth rotating frame and (b) the Earth-Moon rotating frame, with zoomed-in views 
near the Earth appearing in the insets. In Figure 4(a), a circular approximation to the lunar orbit is overlaid 
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Figure 4: Natural motion corresponding to forward propagation of the recent Lunar IceCube deployment 
conditions over 200 days in nondimensional coordinates in (a) a Sun-Earth rotating frame and (b) an Earth- 
Moon rotating frame. 


in black while Figure 4(b) includes the location of the Earth and Moon as gray-filled circles and the libration 
points as red diamonds. As evident in both views represented within this figure, the spacecraft undergoes a 
trailing edge lunar flyby prior to quickly departing the Earth vicinity. Since the Lunar IceCube spacecraft 
must be delivered within a reasonable time interval to its final low altitude lunar orbit to commence scientific 
observations, propulsive intervention is required. 


By introducing low-thrust arcs shortly after deployment, the flyby and post-flyby conditions are sufficiently 
altered to produce a trajectory that remains within the Earth vicinity. To demonstrate the impact of pre- 
perilune low-thrust arcs, consider the following operationally feasible arc sequence: 


e An 8 hr natural arc post-deployment to enable orbit determination, 


e A 2 hr low-thrust-enabled path with the thrust vector aligned with the velocity direction for thruster 
calibration, 


e A 2 hr natural arc for orbit determination and verification of the thruster performance, 


e Three cycles of a 20 hr low-thrust-enabled segment with the thrust vector aligned with the velocity di- 
rection to increase energy and raise the perilune followed by a 10 hr natural arc for orbit determination, 


e Coast until lunar periapsis, and 


e A natural post-flyby arc until that continues until apogee. 


For each low-thrust-enabled arc, a thrust magnitude of 0.9 mN is employed for illustration. This sequence 
of arcs is depicted in Figure 5 in the (a) Sun-Earth rotating frame and (b) Earth-Moon rotating frame. In 
each figure, natural coasting segments are depicted in blue while low-thrust arcs are indicated in green. In 
Figure 5(a), a circular approximation to the lunar orbit is overlaid in black while Figure 5(b) includes the 
location of the Earth and Moon as gray circles and the libration points as red diamonds. The motion of 
the spacecraft is simulated until the first apogee following the lunar flyby. For this example, the apogee is 
prograde with respect to the Earth as viewed in the Sun-Earth rotating frame, 1.e., the motion of the spacecraft 
is instantaneously counterclockwise. In fact, for this particular epoch and deployment state, the incorporation 
of three 20 hr low-thrust arcs with a thrust vector aligned with the velocity direction and separated by coasting 
segments with a duration of 10 hrs raises the lunar periapsis radius to 5080 km from approximately 2314 km 
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Figure 5: Path of Lunar IceCube from deployment to apogee, produced by leveraging both the low-thrust 
engine and a lunar flyby in (a) a Sun-Earth rotating frame and (b) an Earth-Moon rotating frame. 


for the entirely natural trajectory depicted in Figure 4. Note that another deployment state and epoch may 
produce distinctly different behavior.* Nevertheless, in addition to the reduction of speed provided by the 
given finite duration burn, raising the lunar periapsis reduces the acceleration boost supplied by the Moon 
and, therefore, the velocity of the spacecraft after the lunar flyby. As a result, an apogee occurs within the 
Earth vicinity and the path remains, at least temporarily, captured. This apogee is characterized and leveraged 
to link the Earth outbound portion of the transfer to the long phasing and energy adjustment segment. 


By varying the thrust force direction and duration along each of the pre-flyby low-thrust arcs, the location 
and velocity of the first apogee along a temporarily captured path is significantly impacted. Assume the same 
trajectory sequence as in the previous example from Figure 5. While the deployment state is fixed, the di- 
rection of the thrust vector along each of the pre-perilune low-thrust-enabled arcs are varied, impacting the 
apogee conditions. To explore the variation in the state components at this apogee, consider a fixed burn 
duration of 20 hrs. Then, the thrust direction is varied up to 20 degrees relative to the velocity direction and 
is defined as constant as viewed in a Moon-centered VNC frame over each burn segment in a single simu- 
lation. The resulting trajectories are numerically integrated forward in a point mass ephemeris model from 
deployment to the first apogee; the parameters describing both the perilune and apogee are impacted. Figure 
6 depicts the variation in these parameters for various thrust vector directions, where the post-deployment 
lunar flyby conditions are visualized on the b-plane in Figure 6(a) and (b), and the resulting apogees are 
displayed in a Sun-Earth rotating frame in Figure 6(c). In Figure 6(a), each b-plane coordinate describing the 
incoming flyby asymptote is colored by the Jacobi constant of the corresponding perilune state as computed 
in the Earth-Moon CR3BP and indicated via the colorbar. In Figure 6(b), each b-plane coordinate is colored 
by the Jacobi constant of the apogee as computed in the Sun-Earth CR3BP. Then, in Figure 6(c), each apogee 
is colored by the Jacobi constant at apogee, computed in a Sun-Earth CR3BP. Analysis of these figures re- 
veals that pointing the thrust vector up to 20 degrees off the velocity direction can sufficiently reduce the 
energy along the transfer path to ensure temporarily captured motion in the Earth vicinity, while also shifting 
the intersection of the incoming asymptote and the b-plane by thousands of kilometers. Accordingly, the 
perilune radius can also be raised by as much as 2600 km. This variation in the lunar flyby distance signif- 
icantly impacts the radius and energy at the post-perilune apogee, as represented in Figure 6(c). In fact, the 
apogees displayed in Figure 6(c) vary by approximately 760,000 km in distance to the Earth and almost 90 
degrees in their orientation as displayed in a Sun-Earth rotating frame. Furthermore, the Jacobi constant at 
these apogees possesses values between 3.000636 and 3.001185, which bound (above and below) the Jacobi 
constant value of L2, indicating a large variety of the available flow structures in the vicinity of the Earth. 
Accordingly, there may be significant flexibility in the geometry and existence of transfers that return the 
spacecraft to the lunar vicinity. Such insight into the variation in the apogee condition is useful when linking 
the Earth outbound segment to the phasing and energy adjustment segment, as well as the identification of 
feasible transfer regions. 
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Figure 6: Over approximately 729 simulations, impact of varying thrust direction up to 20 degrees from the 
velocity direction on (a) lunar b-plane coordinates colored by the Earth-Moon Jacobi constant at perilune, (b) 
lunar b-plane coordinates colored by the Sun-Earth Jacobi constant at apogee, and (c) post-flyby apogee as 
viewed in the Sun-Earth rotating frame and colored by the Sun-Earth Jacobi constant at apogee. 


Phasing and Energy Adjustment Segment 


Following the post-deployment encounter with the Moon, one option for the Lunar IceCube trajectory 
design is leveraging the gravity of the Sun prior to approaching and capturing into the lunar science orbit. 
While this longest trajectory segment remains within the Earth vicinity, natural dynamical structures within 
the Sun-Earth system can be actively selected to modify both the energy and phasing of the spacecraft in its 
orbit. To reduce the number of deterministic thrusting arcs required along this portion of the Lunar IceCube 
mission trajectory, predominantly natural motions are preferred. Accordingly, Poincaré mapping techniques 
are useful in exploring the geometry of the natural flow corresponding to trajectories that persist within the 
Earth vicinity. These maps are constructed in the CR3BP, a model that offers a reasonable approximation 
to the dynamics in the true Sun-Earth system. Furthermore, these maps allow a prediction for the geometry 
of potentially successful transfers as constructed in this simplified dynamical model, as well as the corre- 
sponding regions of existence for transfers that link to a lunar capture arc in this model and serve as a design 
strategy in the ephemeris model. 


To simplify the visualization of a large array of trajectories at a single energy level in the CR3BP, planar 
apoapsis maps are employed. Construction of these maps within the Sun-Earth CR3BP is demonstrated for a 
sample fixed value of the Jacobi constant, C'; = 3.0008813. At this value of the Jacobi constant, both the L 
and Lz gateways are slightly open and only a small number of trajectories depart the Earth vicinity, enabling a 
clear demonstration of the analysis employed in this investigation. First, feasible initial conditions in the Earth 
vicinity are seeded in the rotating frame using nondimensional configuration space coordinates that lie within 
the zero velocity curves, and between the L; and Lz gateways. These feasible sets of initial conditions are 
also represented as an apoapsis with respect to the Earth. For a state to be considered an apoapsis with respect 
to the Earth in the Sun-Earth rotating frame, it must possess a relative position vector F = [x—1+ py, y, 2], that 
is instantaneously perpendicular to the relative velocity vector, 0 = |, y, 2], such that 7-0 = 0. Furthermore, 
the radial acceleration at an apoapsis must be negative. For this investigation, only planar motion is considered 
when creating apoapsis maps in the CR3BP. Although feasible trajectories within the true ephemeris model 
exist in three dimensional space, the out-of-plane component of motion along each sample transfer is small. 
Accordingly, planar motion in the CR3BP offers a valuable preliminary approximation with the added benefit 
of straightforward visualization. The direction of the velocity vector at each apoapsis is selected across 
the entire map such that the direction of motion is uniformly either prograde with respect to the Earth or 
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retrograde, i.e., counter-clockwise or clockwise, respectively. Thus, for various combinations of the planar 
position components, the direction of the velocity is determined via orthogonality. For a specified value of 
the Jacobi constant, the unit vector along the velocity direction is then scaled using the velocity magnitude, 
computed as v = \/2U* — C. Each initial apoapsis, seeded within the vicinity of the Earth, is then propagated 
forward for a specified number of revolutions about the Earth, from the perspective of the rotating frame. 
Initial conditions that produce trajectories that either impact the Earth or pass through the Sun-Earth L, or 
Lz gateways are discarded. The remaining initial conditions are plotted in configuration space, producing 
a composite representation of the initial apoapses corresponding to trajectories that remain within the Earth 
vicinity as predicted by the Sun-Earth CR3BP. As an example, two apoapsis maps are displayed in Figure 
7 for the Jacobi constant value C'; = 3.0008813, representing trajectories that complete one revolution 
about the Earth without departing through the L; or Lo gateways or impacting the Earth. In Figure 7(a), 
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Figure 7: Apoapsis maps in the planar Sun-Earth CR3BP for trajectories that encircle the Earth once at 
C'; = 3.0008813 for (a) prograde and (b) retrograde initial apoapsis conditions, colored by the angle of the 
final apoapsis after one revolution with the +x direction. 


each initial apoapsis is prograde, while Figure 7(b) displays only retrograde apoapses. For convenience, 
these maps are represented in dimensional Earth-centered rotating coordinates. Gray shaded areas in each 
figure indicate forbidden regions, bound by the zero velocity curves, where motion cannot extend within 
the phase space of the CR3BP for the specified value of the Jacobi constant. Blue points locate apoapses 
that produce trajectories that remain within the Earth vicinity over one revolution and do not impact Earth. 
White regions, however, result in trajectories that do 

not fulfill these criteria. For prograde apses, these y 

white regions are bound by the stable manifolds of 
the L; and Lz Lyapunov orbits, while the white re- 
gion in Figure 7(b) corresponds to trajectories that yee yee 
pass inside the radius of the Earth.!9 Furthermore, — To Sun 
red diamonds locate the equilibrium points, while a §=<— 
light blue circle indicates the location of the Earth 

and the black circle represents a circular approxi- x-l+u<0 
mation to the orbit of the Moon. Then, each initial y<0 
apoapsis on these maps is colored by the angle of 
the final apoapsis after one revolution in the rotating 
frame relative to the Earth-L line. A legend for this 
color scheme appears in Figure 8, with the Earth lo- 
cated at the center as a blue-filled circle. To mitigate 
the effect of complicating the map visualization, the colors in this scheme vary every 30 degrees, indicating 
whether the final apse along a trajectory arc lies in the center of a quadrant, or close to a neighboring quadrant. 
For instance, initial conditions on an apoapsis map are colored teal, blue or dark blue if the final apse along 


x-1l+u<0 x-l+u>0 
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Figure 8: Color scheme used in the mapping 
strategy, indicating the angular location of the final 
apogee relative to the Sun-Lo line. 
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the trajectory, propagated only for a specified number of revolutions around the Earth, possesses a state in the 
Sun-Earth rotating frame such that x -—1+ p< Oand y < 0, 1.-e., it falls within the bottom left quadrant of the 
Earth-centered view. While this color scheme for the initial conditions on an apogee map does not capture 
the exact location of the final apses along the associated paths, it does enable rapid and straightforward eval- 
uation of candidate transfers to identify those to be analyzed in further detail. The insight obtained through 
this application of Poincaré mapping strategies may guide the numerical targeting of outgoing lunar flyby 
conditions that subsequently deliver the Lunar IceCube spacecraft on a natural transfer path that requires lit- 
tle propulsive effort. For instance, white regions should not be targeted to avoid departing the Earth vicinity 
or impacting the Earth, while certain colored regions and their size may indicate the geometry of feasible 
transfers and their sensitivity to uncertainty in the initial conditions. In a chaotic and nonlinear system, sim- 
ply integrating the finite set of post-deployment perilune conditions in Figure 6 forward for multiple Earth 
revolutions would yield only limited results that do not reflect the range of feasible transfer geometries, or 
supply insight into the corresponding regions of existence. However, mapping strategies supply a capability 
to visualize and rapidly identify natural candidate arcs in a simplified, yet representative, dynamical model. 
Such insight enables a well-informed process for constructing an initial guess for a trajectory that delivers 
the spacecraft from a highly energetic deployment state near the Earth to the vicinity of the Moon, as well as 
a rapid exploration of the trajectory design space. 


Lunar Capture Segment 


To locate arcs that approach the Moon and eventually result in captured motion, techniques from the 
construction of invariant manifolds are leveraged.* Recall that the Lunar IceCube spacecraft is required to 
collect observations from a low altitude lunar orbit that is both elliptical and highly inclined. One strategy for 
reaching this orbit from a state that is located well beyond the lunar radius is based on the design of a trajectory 
that passes through the Ly gateway with a Jacobi constant that is below the value corresponding to the Lo 
libration point, i.e., C'y(L2). Then, the low-thrust engine is activated in a direction close to the anti-velocity 
vector, effectively reducing the spacecraft velocity. Simultaneously, with an active thrust force, the Jacobi 
constant increases, with the Lz and, then, L; gateways closing, thereby bounding the motion of the spacecraft 
to the vicinity of the Moon and yielding captured motion. By modifying the thrust direction over a long time 
interval, the final science orbit is accessed. A significant challenge in designing a trajectory to deliver the 
spacecraft from a highly energetic deployment state to the final science orbit via this strategy is locating a 
trajectory that passes through the Lz gateway. Since the Sun-Earth-Moon system is inherently nonlinear 
and chaotic, an uninformed search for trajectories that achieve this goal may be cumbersome, tedious and 
unsuccessful. However, techniques from manifold computation approaches are useful: a desired lunar orbit 
is discretized and each fixed point is numerically integrated backwards in time in a point mass ephemeris 
model, with the thrust vector directed in the anti-velocity direction.* In this investigation, to construct a 
successful lunar capture, one of the significant hurdles in trajectory design for the Lunar IceCube spacecraft 
is identifying trajectories that deliver the spacecraft to the lunar vicinity with a Jacobi constant value above 
that of L,; and Lz. To achieve this goal, it is necessary to construct an initial guess for a lunar capture arc 
that completes multiple nearly-polar revolutions around the Moon. Higher order gravitational models, which 
significantly influence motion close to the bodies, can be then be employed in additional analyses to exactly 
target motion that subsequently reaches the low lunar orbit. 


To explore the dynamical flow toward the destination orbit, a desirable lunar orbit is defined and integrated 
in negative time with the low-thrust engine activated. The desired science orbit, with a Jacobi constant above 
that of Ly and Lo, is described by a semi-major axis of 4,287 km, an eccentricity of 0.5714, an inclination 
of approximately 90 degrees, and a predefined final epoch. The right ascension of this orbit is sampled 
within the range 2 = [0,360). States with true anomalies in the range 6* = [0,360), are sampled along 
this orbit and integrated backwards in time in a point mass ephemeris model with the thruster activated until 
reaching the Jacobi constant of £;. With the thruster still activated, the trajectory is integrated further until 
the subsequent apse (periapsis or apoapsis). Departing from this state, the thruster is off and the spacecraft is 
naturally propagated further backwards in time until it reaches the L; or Lz gateways, impacts the Earth or 
until reaching a maximum integration time. The parameters describing this trajectory are stored if they pass 
through the Lz gateway and encounter an apogee beyond the lunar radius within a reasonable time frame. 
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This process is then repeated by integrating backwards in time with the low-thrust engine activated to the 
next apse until the trajectory passes below the Earth radius or departs the Earth vicinity. A sample trajectory 
sequence is displayed in Figure 9 for the prescribed values of a, e, and 7, a final epoch of 28601 MJD, a 
final mass of 13.5 kg, a RAAN of 150 degrees and a true anomaly of 135 degrees. This sample trajectory 
is visualized in (a) a Moon-centered inertial J2000 equatorial frame and (b) a Moon-centered Earth-Moon 
rotating frame. In both plots, green indicates thrusting segments while blue curves identify natural coasting 
arcs. To generate a set of trajectories that may be used to construct a viable initial guess, a final spacecraft 
mass at arrival to the science orbit must be set. Since this quantity is not known a priori, it is estimated 
iteratively using the average thrust time for the spacecraft to depart the L2 gateway. This process, 1.e., 
generating a set of sample capture arcs, is demonstrated for the selected science orbit and a final epoch of 
28695.0 MJD, with values for 2 and @* each varied in increments of 45 degrees, respectively. The final 
spacecraft mass is estimated to be equal to 13.5 kg. Trajectories that experience an apogee within 250 days 
of passing through the Ly gateway are stored and plotted in Figure 10(a) in the Sun-Earth rotating frame. 
Then, the recorded apogee for each of the stored trajectories is indicated by a circle colored by the epoch 
in terms of modified Julian date in Figure 10(a). These trajectories exhibit a variety of characteristics and 
possess apogees that occur in various regions of the phase space. To aid visualization and comparison to the 
phasing and energy adjustment segment, these apogees are also depicted on an apoapsis map in Figure 10(b). 


(a) Moon-Centered Inertial J2000 Frame (b) Moon-Centered Earth-Moon Rotating Frame 
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Figure 9: Sample arc that captures into the desired lunar science orbit, displayed in (a) a Moon-centered 
inertial J2000 equatorial frame, and (b) a Moon-centered Earth-Moon rotating frame. Green arcs indicate 
that the thruster 1s activated in the anti-velocity direction, while blue arcs correspond to natural motion 


(a) Epoch (MJD) (5) Epoch (MJD) 
x 10° aie x 10 x10. 
Earth 
2.852 1 « le 2.852 
oe ee . 
a 2.85 a ae 2.85 
> > ame 4 
2.848 | 2.848 
-10 -5 O 5 -10 -5 O 5 
x (km) y 19° ea 0) 


Figure 10: (a) Trajectories that deliver the spacecraft to the lunar vicinity as viewed in configuration space in 
the Sun-Earth rotating frame and (b) representations of these paths on an apoapsis map colored by the apogee 
epoch in MJD. 
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Each apogee, represented by a planar projection of the configuration space variables in the Sun-Earth rotating 
frame, is colored by its corresponding epoch to identify lunar approach arcs that originate at a time similar to 
the time associated with the final state along the previous segment. Of course, this map does not represent all 
parameters describing the apogees within this set. However, this reduced dimension mapping does enable a 
rapid and straightforward evaluation of the trade space. From this set of apogees, candidate paths that deliver 
the spacecraft to the lunar vicinity are straightforwardly identified and incorporated into the construction of 
an initial guess for a complete end-to-end trajectory design. 


Sample Trajectory Construction 


Arcs comprising each segment of the transfer sequence are selected using apoapsis maps to construct an 
initial guess for a path that delivers the spacecraft from an initial deployment state to the lunar vicinity and is 
then corrected in a low-thrust-enabled point mass ephemeris model. This process is demonstrated for a single 
transfer using the most recent deployment information (current as of December 2016). As the deployment 
state and epoch continues to evolve throughout the EM-1 mission planning phase, the transfer sequence and 
Poincaré map definition may require updates; however, the general framework, constructed and demonstrated 
in this investigation, remains useful. 


An Earth outbound transfer arc influences the potential geome- 


tries for feasible transfers in the phasing and energy adjustment seg- 5 

ment. For this example, a single sequence of coasting and thrusting g x 10 

arcs is selected using a cyan point from Figure 6. Consider a sample a 

trajectory generated in a point mass ephemeris model of the Earth, 6} 

Moon and Sun and plotted as a planar projection in Sun-Earth ro- Lunar 

tating coordinates in Figure 11. Recall that blue arcs identify nat- 4) Orbit 

ural motion while green arcs indicate that the 0.9 mN low-thrust “@ 

engine is activated and a black circle approximates the lunar or- ~~ 2 

bit. The final apogee along this transfer segment is prograde when e 

viewed in the rotating frame, possessing a Jacobi constant value of 0 

approximately C'; = 3.000854 and an epoch of 28436.5 MJD. Fur- 9} 

thermore, at the end of this segment, the spacecraft mass is equal 

to 13.99 kg. Each of these characteristics significantly impact the _At 

construction of the remainder of the trajectory. 6 4 2 0 2 4 
To identify a feasible path for the energy and phasing adjust- x (km) X 10 


ment segment, the Poincaré mapping strategy is employed. For 
this example, consider transfer arcs that encircle the Earth once, es- 
sentially reducing the time of flight along this portion of the path. 
An apoapsis map is then constructed for this type of transfer at 
C'; = 3.000854 with prograde initial conditions, matching the Ja- 
cobi constant value and the direction of motion at the apogee that 
designates the end of the Earth outbound segment. The resulting 
map is displayed in Figure 12. Recall that initial apogees that produce planar trajectories encircling the Earth 
once without departing through the L,; or Ly gateways or impacting the Earth are plotted on this map in 
Earth-centered rotating coordinates. Each initial condition is colored, according to the legend displayed on 
the right of Figure 12, by its fate after a subsequent revolution, represented by the angle of the final apoapsis 
state along the trajectory relative to the Sun-L line. In addition, a circular approximation to the lunar orbit 
is plotted in black and the Earth is located by a gray-filled circle; the equilibrium points are plotted as red 
diamonds. Then, forbidden regions are again indicated by gray shaded regions. Overlaid on this apoapsis 
map is the location of the apogee state at the end of the Earth outbound segment, plotted as a gray star. Recall 
that the Earth outbound transfer arc is computed in a three-dimensional point mass ephemeris model, while 
the trajectories in this energy and phasing adjustment segment are propagated in a planar Sun-Earth CR3BP 
model. Thus, an exact match is not expected between the initial condition associated with the actual transfer 
arc in this segment and the end of the Earth outbound segment. However, the CR3BP does supply a reason- 


Figure 11: Planar projection of the 

selected outbound segment in Earth- 
centered rotating coordinates in the 
Sun-Earth frame. 
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Figure 12: Apoapsis map in the planar Sun-Earth CR3BP for trajectories that encircle the Earth once at 
C'; = 3.000854 for prograde initial conditions, colored by the angular location of the final apoapsis as 
indicated in the legend. The final apogee along the Earth outbound segment is indicated by a gray star. 


able approximation to the true dynamics in the Sun-Earth system. Selecting an initial apogee near the final 
apogee along the previous segment should mitigate any necessity for a long thrusting arc to connect the two 
segments, while still providing sufficient flexibility to select the desired transfer geometry. Thus, this initial 
condition could lie in either a yellow, green or orange region, each of which yield transfers that possess a 
final apogee near the +y-axis. In this example, the location of the final apoapsis along this transfer arc offers 
a rapid assessment of the approximate geometry for the transfer. This coloring scheme is also useful if the 
epoch at which the spacecraft must reach the final science orbit is constrained. An initial condition in the 
yellow region near the gray star is selected close to the boundary of the orange area and propagated in the 
CR3BP to produce a planar trajectory that is plotted in the Sun-Earth rotating frame in Figure 13(a). This 
transfer, with a time of flight equal to 83.3 days, is then discretized into 20 subarcs, equally spaced in time. 
The spacecraft state at the beginning of each arc, as computed in the CR3BP, is then propagated forward in 
the point mass ephemeris model of the Sun, Earth and Moon. In support of this simulation, an initial epoch, 
equal to the epoch at the end of the previous segment, 28436.5 MJD, is assumed. As an estimate, the final 
apogee at the end of this transfer arc is approximated as 83.3 days following the initial epoch, 1.e., 28519.7 
MJD. Note that, since only natural motion is leveraged along this segment, the spacecraft mass remains un- 
changed at the value of 13.99 kg. The resulting sequence of arcs is plotted in planar configuration space in 
Figure 13(b) with the nodes at the beginning of each segment located with green circles. Of course, each 
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Figure 13: Selected transfer arc for energy and phasing adjustment segment plotted in the Sun-Earth rotating 
frame: (a) propagated in the CR3BP and (b) discretized and integrated in the point mass ephemeris model. 
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of these arcs is three-dimensional and no longer produces a single continuous path. However, these nodes 
serve as an initial guess for the Sun-Earth segment along the complete transfer, one that is later adjusted in 
a corrections algorithm. Additionally, the map in Figure 13 supplies insight into the available geometries of 
candidate arcs for the Sun-Earth segment. For instance, the thin colored regions in the vicinity of the gray 
star in Figure 12 indicate that the location of the final apogee along the Sun-Earth segment of the transfer may 
be sensitive to the initial apogee state: such sensitivity may be beneficial in navigating back to a reference 
science orbit approach path with minimal propellant usage following small uncertainties in the deployment 
information, thruster performance or flyby conditions. However, a sufficiently large uncertainty that places 
the initial apogee close to the small white region in Figure 12, embedded within the light blue structure in the 
top left of the figure, may necessitate the use of long thrust arcs to ensure that the trajectory does not depart 
the Earth vicinity or impact the Earth. 


A final transfer arc is, concurrently, selected to deliver the spacecraft to the lunar vicinity. Recall that the 
transfer arc selected in the energy and phasing adjustment segment produces a final apogee that is shifted to 
the Lz side of the Earth in a Sun-Earth rotating frame. Thus, lunar approach arcs that possess an initial apogee 
on the [2 side of the Earth in the Sun-Earth system with an epoch close to 28519.7 MJD are sought. Although 
a close match in energy between the beginning of the lunar approach segment and end of the phasing and 
energy adjustment segment is sought, the value of C’'; is no longer constant in the ephemeris model, especially 
with the incorporation of an additional gravitational contribution. Thus, a prediction of the value of C’; at 
the end of the energy and phasing adjustment segment from the CR3BP may not be truly reflective of the 
actual value when corrected in an ephemeris model that also incorporates lunar gravity. Accordingly, the 
Jacobi constant is only loosely employed as a filter for selecting the apogee that defines the lunar approach 
arc. Alternatively, the epoch, location and direction of motion are incorporated as the dominant selection 
parameters. To produce an apoapsis map for identifying a lunar approach arc, the process outlined earlier is 
employed along with the specified science orbit. An appropriate final epoch at the science orbit is iteratively 
determined and set equal to 28693 MJD, while the right ascension for the orbit and the true anomaly for 
discretization are sampled within the range [0, 360) in increments of 40 and 10 degrees, respectively. The 
final spacecraft mass is set to 13.5 kg. The resulting trajectories, in particular, those that pass through the Lo 
gateway in backwards time, are stored and plotted in Figure 14(a) and (b) in the Earth-centered Sun-Earth 
rotating coordinates. In Figure 14(a), each apogee is colored by the number of days past the desired epoch, 
1.e., the epoch at the end of the previous transfer segment and equal to 28519.7 MJD, while in Figure 14(b), 
each apogee is colored by the Jacobi constant as computed in a Sun-Earth CR3BP. Overlaid on these plots is a 
gray star which locates the final apogee of the previous arc. Since arcs selected from this map should possess a 
position, epoch and energy similar to those of the reference apogee, a lunar approach apogee close to the gray 
star and in the red region of Figure 14(a) is preferred. This apogee should also possess a Jacobi constant close 
to the value corresponding to the final apogee along the Sun-Earth segment, i.e. C = 3.000854, thereby also 
lying in a light blue to green region of Figure 14(b). If no apogees are detected in the desired regions of the 
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Figure 14: Apoapsis map representing trajectories that approach the lunar science orbit, projected into planar 
configuration space and colored by the epoch at each apogee. 
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apoapsis maps near the final apogee along the previous segment, the final epoch can be iteratively modified 
and a new map constructed. Of course, the difference between the parameters describing the apogee at the 
end of the phasing and energy adjustment segment and the apogee at the beginning of the lunar approach arc 
impact the quality of the initial guess and, therefore, the potential for the differential corrections process to 
identify a continuous solution. Fortunately, in the current example, candidate apogees lie within desirable 
areas near the gray star. Thus, this map is employed to define a small range of RAAN and true anomaly 
values and iteratively recomputed to identify a good candidate arc: one corresponding to a final science orbit 
epoch of 28693.01 MJD, RAAN = 199.2 deg, TA = 68 deg. A lunar approach arc is selected in this vicinity, 
producing the trajectory plotted in Figure 15 and viewed in: (a) an Earth-centered view in the Sun-Earth 
rotating frame, (b) a planar zoomed-in view along the final portion of the transfer in the Moon-centered 
Earth-Moon rotating frame, and (c) a three-dimensional view of the final portion of the transfer in the Moon- 
centered Earth-Moon rotating frame. Along the trajectory in Figure 15, green segments indicate low-thrust 
arcs, while a blue line identifies a natural arc. Using this figure, the selected initial condition occurs at an 
epoch equal to 28523.6 MJD with a Jacobi constant value C'; = 3.000917 and a spacecraft mass equal to 
13.92 kg; this mass is approximately 0.07kg below the spacecraft mass at the end of the previous transfer 
segment and the epoch is discontinuous by approximately 4 days. In addition, the incoming lunar capture arc 
exhibits a lunar flyby distance of 100 km prior to activating the low-thrust engine. Of course, these parameters 
are adjusted during the iterations in a corrections algorithm, with short low thrust arcs introduced to assist in 
bridging the discontinuities. Since the goal of this investigation is the construction of a transfer that delivers 
the spacecraft to the lunar vicinity, the selected trajectory is clipped to include only a few revolutions of the 
Moon, rather than the entire lunar approach arc. In fact, these revolutions around the Moon serve to guide 
the corrector to a trajectory that reaches the lunar vicinity without the extra computational time associated 
with correcting a longer path. The exact trajectory to the final lunar science orbit is more appropriately 
identified in subsequent analyses that incorporate higher fidelity models to also capture perturbations from 
lunar gravitational harmonics. 


Each of the arcs selected via dynamical systems techniques enables the construction of an initial guess 
for a trajectory that is ultimately corrected in a low-thrust-enabled point mass ephemeris model. In par- 
ticular, each of the three segments depicted in Figures 11, 13 and 15 are concatenated chronologically 
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Figure 15: Selected trajectory that approaches the lunar vicinity: (a) a view of the entire trajectory in an 
Earth-centered Sun-Earth rotating frame, (b) zoomed-in planar view of the final portion in a Moon-centered 
Earth-Moon rotating frame and (c) a three-dimensional view in the Earth-Moon rotating frame. 
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with short low thrust arcs introduced between each segment, 1.e., the Earth outbound segment is followed 
by a short low thrust arc to connect to the energy and phasing adjustment segment, and then an addi- 
tional low thrust path is followed by the lunar approach arc. Recall that the corrections algorithm em- 
ployed in this investigation is formulated for simultaneous adjustment of a sequence of nodes described 
in terms of the body-centered position and velocity coordinates in a J2000 inertial frame, nondimensional- 
ized by the characteristic quantities corresponding to the Sun-Earth or Earth-Moon system, as appropriate. 
Additional information that is recorded at each node includes 

the nondimensional time past a reference epoch (e.g., the de- ale Discontinuities in 
ployment epoch), the integration time along the arc of inter- 10! a siehe. (ime cnass 
est as well as the spacecraft mass. For low-thrust-enabled 

arcs, the three components of the thrust direction unit vector 
are also stored as variables. To convert the initial guess in 
Figure 16 into a sequence of nodes, natural and low-thrust- 
enabled segments are first separated. Then, each segment 
is discretized by a user-defined number of nodes that are 
equally spaced in time. For this sample trajectory, each of 
the states along the discretized initial guess is defined relative 
to the Moon. In addition, three constant vectors supply the 5 0 5 

corrections algorithm with information that is vital to the nu- —- Lunar Flyby x (km) eT 

merical integration process and evaluating continuity across 

all arcs, as well as accurately incorporating the low-thrust Figure 16: Discontinuous initial guess for a 
engine when appropriate. Specifically, these vectors reflect complete trajectory, plotted in Earth-centered 
activation of the low-thrust engine, the basepoint body used Sun-Earth rotating coordinates. 

for defining the VNC unit vectors describing the thrust direc- 

tion, and the body used as a basepoint for defining the state 

information. This initial guess is displayed in Earth-centered Sun-Earth rotating coordinates in Figure 16. 
Recall that green arcs locate low-thrust-enabled segments while blue indicates natural motion. Although the 
free variables representing each state are formulated in the inertial frame when inputted to the corrections 
algorithm, the Sun-Earth rotating view offers an intuitive representation for verification of the construction of 
the initial guess. Note that this initial guess exhibits clear discontinuities in the state, time and mass variables 
between each phase of the trajectory. However, short low-thrust arcs are added at these discontinuities to 
support the identification and computation of a continuous path. While selecting different arcs may reduce 
these discontinuities, this initial guess serves as an interesting example for demonstration purposes. 


y (km) 


The constructed initial guess is supplied to a corrections algorithm that leverages a low-thrust-enabled 
point mass ephemeris model. Since this complete trajectory is comprised of arcs that pass close to each of 
the primaries as well as those that exist beyond the lunar radius, numerical issues due to increased sensitivity 
and the presence of significant discontinuities may appear. To mitigate such problems, a modified Newton’s 
method is employed to update the free variable vector. Recall that the update equation for a Newton’s method 
is formulated under the assumption that the current guess for a variable is close to the true solution. However, 
for an initial state that is not sufficiently close to a true solution for the constraint vector, the Newton strategy 
may diverge. Thus, an attenuation factor, y, is introduced to the update equation for an underdetermined 
system such that, at each step, the free variable vector is iteratively modified as follows: 


Vil — 7 _ yDF(V)? DEV). DEV) 1 F(V) ©) 


for a free variable vector at the i-th iteration, V*, a constraint vector, F’, and a rectangular Jacobian matrix, 
DF. Note that 7 can be assigned any value between 0 and 1. Once the magnitude of the constraint vector 
falls below a threshold value, the update equation reverts to the traditional Newton’s method, 1.e., y = 1. 
The exact value for both + and the threshold on the constraint vector magnitude may be tuned by the user as 
the performance of the algorithm is evaluated. Alternative methods to eliminate any issues associated with 
ill-conditioning would serve as an interesting avenue of investigation in future analyses. 


Corrections on the initial guess in Figure 16 enable the detection of a continuous trajectory in a low-thrust- 
enabled ephemeris model that follows the constructed itinerary and delivers the spacecraft from a highly 
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energetic deployment state to the lunar vicinity. The state, time and mass continuity constraints, as well as 
the positive integration time and unit thrust vector direction constraints, are simultaneously applied in a cor- 
rections algorithm with updates implemented until the trajectory is continuous to within a nondimensional 
tolerance on the constraint vector equal to 10~1°. Following the completion of this iterative process, a cor- 
rected trajectory is identified. Additional constraints on the initial condition along the path are also enforced 
to match the specified deployment conditions. This solution is displayed in Figure 17 in (a) Earth-centered 
Sun-Earth rotating coordinates, (b) a zoomed-in view of the lunar vicinity in Moon-centered Earth-Moon 
rotating coordinates, (c) Earth-centered J2000 equatorial inertial coordinates and (d) Moon-centered J2000 
equatorial inertial coordinates. Recall that green lines indicate low-thrust-enabled arcs while blue curves 
indicate natural motion. Comparison to the initial guess as constructed in Figure 16 reveals that the gen- 
eral geometry of the transfer is retained, demonstrating the value of leveraging natural dynamical structures 
within multi-body systems in the construction of a complex path for a low-thrust-enabled CubeSat. Although 
the final portion of the path does not capture the entire lunar approach path, the goal of this investigation is 
certainly achieved: recovery of a path that delivers the spacecraft from a highly energetic deployment state to 
the lunar vicinity, while maintaining a geometry consistent with the constructed initial guess. Identification 
of a trajectory that eventually encircles the Moon, with a sufficiently high inclination and an energy at which 
the L; and Lz gateways are closed, addresses the most significant challenges in the construction of an initial 
guess for a feasible end-to-end transfer. Additional analysis in a higher-fidelity model that more accurately 
reflects the lunar gravitational field can be leveraged to identify subsequent low-thrust-enabled arcs that reach 
the desired lunar science orbit. Furthermore, corrections and optimization processes may be applied to this 
solution to reduce the total thrusting time and, therefore, propellant usage. This process may be used in 
future analyses to rapidly produce a variety of transfers, supporting an exploration of the available transfer 
geometries and corresponding regions of existence.!° 
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Figure 17: Continuous sample trajectory (a) plotted in Earth-centered Sun-Earth rotating coordinates, (b) a 
zoomed-in view in Moon-centered Earth-Moon rotating coordinates, (c) Earth-centered J2000 inertial coor- 
dinates and (d) Moon-centered J2000 inertial coordinates. 


CONCLUDING REMARKS 


The insight gained throughout this investigation demonstrates that a trajectory design framework, con- 
structed using dynamical systems techniques, is valuable in designing a complex path for a spacecraft with 
limited propulsive capability within the chaotic gravitational environment of the Earth, Moon and Sun. The 
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continuous solution constructed here serves as a suitable initial guess for corrections to a complete end-to-end 
path in an operational-level modeling environment and, potentially, an optimization algorithm. Of course, the 
exact transfer sequence and map definition may require updates as the deployment state frequently evolves 
throughout development of the mission. However, the framework developed in this investigation remains 
useful for constructing trajectories for a small spacecraft that must be delivered from a high energy deploy- 
ment state near the Earth to a low lunar orbit. In fact, the design framework enables rapid redesign of this 
trajectory as the deployment information or spacecraft hardware specifications evolve by mitigating the chal- 
lenges associated with searching for a point solution in a chaotic dynamical regime, while also supporting the 
exploration of a large and complex design space. 
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